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Introduction 


During  the  recent  decade,  iron  has  been  the  object  of  extensive  study  by 
various  first  principles  methods  (for  references  and  comparison  of  results  for 
iron  obtained  by  different  methods  see,  e.g.  [1-4]).  However,  the  self-  consistent 
spin-polarized  calculations  on  BCC  iron  inevitably  failed  to  predict  the  relative 
stability  of  the  ferromagnetic  (FM)  BCC  phase  with  respect  to  the  nonmagnetic 
(NM)  FCC  (see  discussion  in  [4]).This  failure  is  believed  to  be  mainly  due  to  the 
local  spin-density  approximation  used  in  the  calculations. 

Because  of  this  fundamental  difficulty,  at  the  present  time,  an  entirely  ab 
initio  analysis  of  structural  phase  transformations  in  iron  is  impossible.  However 
calculations,  which  do  not  require  comparison  of  energetics  of  different  phases, 
like  modeling  of  grain  boundaries  and  free  surfaces,  can  be  successfully 
performed  [5-7].  An  alternative  approach  to  the  totally  ab  initio  analysis  is  to 
introduce  into  the  theory  an  adjustable  parameter  in  order  to  make  calculations 
more  consistent  with  experimental  observations. 

Along  this  path,  in  Refs.  [8-10],  we  have  chosen  to  calculate  the 
equilibrium  magnetic  moments  as  well  as  the  magnetic  contributions  to  the 
ground  state  energies  of  iron  using  the  Stoner  model  of  itinerant  ferromagnetism 
[11],  rather  than  performing  spin-polarized  calculations.  The  Stoner  exchange 
parameter,  /,  can  then  serve  as  an  adjustable  parameter.  Having  made  only  one 
adjustment ,  such  a  procedure  enabled  us  to  perform  the  complete  analysis  of  the 
relative  stability  of  BCC  and  FCC  phases,  as  well  as  the  energetics  of  the  BCC- 
FCC  lattice  deformation  in  iron. 

In  spite  of  the  tremendous  success  in  development,  in  the  recent  decade,  of 
new  first-principles  methods  of  modeling,  as  well  as  efficient  computer  codes  and 
proliferation  of  powerful  high-speed  computers,  a  totally  first-principles 
modeling  of  important  processes  fundamentally  affecting  the  mechanical 
properties  of  metallic  alloys,  is  too  computationally  intensive  and  therefore 
practically  impossible. 

In  response  to  the  pressing  necessity  of  filling  this  void,  the  efficient  semi- 
empirical  methods,  such  as  different  versions  of  the  so-called  Embedded  Atom 
Method  (EAM)[12],  or  Angular  Force  Method  (AFM)[13],  have  been  developed. 


However,  the  existing  semi-empirical  methods  for  use  in  atomic  computer 
simulations  in  iron  still  have  not  explicitly  taken  account  of  the  ferromagnetism 
per  se.  Rather,  the  adjustable  parameters  just  have  been  fitted  to  mechanical  and 
thermal  properties  of  iron  At  the  same  time,  when  the  magnetic  properties  are 
essential  in  affecting  the  atomic  behavior  (e.g.  in  the  vicinity  of  crystal  defects), 
these  methods  may  be  inadequate. 

To  remedy  the  situation,  we  suggest  the  procedure  of  explicitly  calculating 
the  magnetic  contribution  to  the  energy  of  individual  atoms,  as  a  function  of  their 
atomic  environment.  This  procedure  is  based  on  the  Stoner  theory  of  itinerant 
ferromagnetism,  and  uses  the  information  that  has  been  accumulated  in  our 
earlier  research  on  iron,  mentioned  above. 

The  plan  of  this  paper  is  as  follows.  In  Sec.  2  we  discuss  the  Stoner  model 
of  itinerant  ferromagnetism.  In  Sec.  3  ,we  discuss  the  parametrization  of  the  the 
ingredient  quantities  to  be  used  in  semi-empirical  calculations.  Sec.  4  summarizes 
the  procedure  and  discusses  potential  applications  of  the  new  method  to  EAM  and 
AFM. 


The  Stoner  Model  of  Itinerant  Ferromagnetism 

The  Stoner  theory,  first  suggested  in  1939  [11],  has  been  successfully  used 
in  recent  years  in  estimating  both  the  equilibrium  magnetization  and  magnetic 
energy  of  band  electrons.  This  was  made  possible  as  a  result  of  a  rigorous 
formulation  of  the  Stoner  model  as  a  perturbation  approach  in  terms  of 
microscopic  electronic  theory  [14-17].  Particularly,  the  fundamental  parameter  of 
the  theory,  the  Stoner  exchange  parameter,  /,  was  understood  in  terms  of 
density-functional  characteristics.  In  iron,  the  theory  explained  the  metamagnetic 
behavior  of  the  FCC  phase  [14,18,7].  The  Stoner  approach  in  combination  with 
self-consistent  non-spin-polarized  calculations  enables  one  to  perform  the  detailed 
analysis  of  ferromagnetic  (FM)  behavior,  as  well  as  identify  all  the  possible 
magnetic  stationary  phases;  both  stable,  metastable  and  even  unstable,  and  find  the 
areas  of  their  emergence.  Such  an  analysis,  using  traditional  spin-polarized 
calculations  is  at  present  either  too  cumbersome  and  practically  impossible. 

The  Stoner  model  in  its  original  formulation  postulates  that  the  change  of 
energy  upon  forming  a  FM  state  with  moment  m  consists  of  two  parts.  The  first 
is  the  exchange  energy  contribution,  -1/4  /  m^,  where  the  exchange  parameter  7 
is  a  constant.  The  second,  the  kinetic  energy  term  is  found  by  forming  two 
subbands  for  spin  up  and  down  electrons  by  flipping  m/2  spin  down  electrons 
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from  just  below  the  non-magnetic  (NM)  Fermi  level  into  the  unoccupied  spin-up 
states  just  above  the  Fermi  level.  As  was  shown  in  [14-17],  this  procedure 
corresponds  to  the  first  order  perturbation  theory  in  m/ny  (  ny  is  the  number  of 

valence  electrons  per  atom).  Thus,  for  a  given  m,  the  magnetic  contribution  to 
the  total  energy  is: 

m 

Em=  1/2  j  m'  /N(m')  dm'  -  1/4  I  m2  (1) 

0 

where  N(m)  is  the  NM  density  of  states  averaged  between  the  Fermi  levels  of 
spins  up  and  down  electrons  as  found  from  the  rigid  subband  shift.  The  procedure 
of  'constructing'  N(m)  is  illustrated  in  Fig.  1. 


NM  d-DOS 


Fig.  1.  Schematic  construction  of  the  N(m)  function. 

The  stationary  state  requirement  : 

d  Ejn  /3m  =0, 

gives,  apart  from  the  "trivial"  solution,  m=0,  the  criterion  for  arising  a  FM  state: 

I  N(m)  =  1  (2) 

Suppose  Eq.  (2)  has  a  solution,  m.  Then  the  corresponding  FM  state  is 
stable  (  32Ein  /3m  2>  0)  if 

3N(m)/3m<0  (3) 
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otherwise  it  is  unstable  (  /dm  2<  0).  However,  even  if  Eqs.  (2,3)  hold,  the 
EM  state  may  not  occur  if  E^  >  0.  In  this  case  the  EM  state  is  metastable. 

As  for  the  NM  state,  m=0,  it  is  stable  and  may  coexist  with  a  EM  state 
(metamagnetic  situation),  only  so  far  as 

[  32Em  /3m  2  >  o  (4) 

or,  equivalently, 

I  N(Ef)  <  1  (5) 

where  N(Ef)  is  the  NM  density  of  states  (  DOS)  at  the  Eermi  level  (N(Ef)=N(0)). 

Erom  the  perturbation  theory  analysis  [14-17]  ,  /  can  be  found  in  terms  of 
the  NM  system.  Erom  the  linear  response  theory  [17]  it  follows; 

/  =  jd3ry2(r)|K(r)l  (6) 

(the  tree-dimentional  intergal  is  taken  over  the  whole  volume  of  crystal)  where 


and 


Y(r) 


Ip-E.)  I\|/(r)l^ 

N(Ep) 


K(r)=  1/2  [  d  2  Exc(r,m)  /  dm  2  ] 


Here  Ei,  \|/i  (r)  are  respectively  the  eigenvalues  and  the  wave  functions  of  the  NM 
system;  Exc(r7iii)  is  the  exchange-correlation  functional. 

An  important  property  of  of  the  averaged  DOS  N(m)  is  that,  being 
multiplied  by  the  d-electron  "band  width",  W  ,  it  happens  to  be  virtually 
independent  of  volume,  or  the  Wigner-Seitz  (WS)  radius,  s.  Then  Eq.  (2)  reads: 


A(m)=W(s)//(s) 


(7) 


Where  A(m)  =  W*N(m).  Eig.  2  shows  N(m)  curves  for  two  s  values.  The  right- 
hand  side  of  the  above  Eq.  does  depend  on  s,  and  this  dependence  is  important. 

Until  recently,  the  Stoner  parameter,  I ,  has  been  believed  to  be  essentially 
a  constant,  independent  of  both  the  volume  and  the  crystal  structure  of  the  metal. 
Our  calculations  on  both  the  BCC  and  ECC  iron  [8-10]  revealed  the  monotonic. 
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Fig.  2.  A/^(m)s  for  BCC  lattice  for  two  WS  radii. 

though  rather  weak,  dependence  of  I  on  the  WS  radius,  s.  Both  the  magnetic 
energy,  Eq.(l),  the  equilibrium  magnetic  moment,  m,  (as  found  from  Eq.  (2)), 
and  the  equilibrium  atomic  volume,  happen  to  be  rather  sensitive  to  the  values  of 
I.  The  idea  of  papers  [8-10]  was  to  adjust  the  value  of  /  ,  so  that  the  equilibrium 
WS  radius,  Sq,  for  the  FM  BCC  phase  be  equal  to  the  experimental  value.  In 
Table  I  we  compare  the  calculated  Sq,  for  /  =(3  Iq  [  where  Iq  is  the  "ab  initio” 

value,  Eq.  (6),  and  p=  1.000, 1.025,  1.050,  1.075  and  1.090].  One  can  see  that  for 
P=  1.075,  the  equilibrium  WS  radius,  So=2.659  a.u.,  almost  matches  the 

experimental  value.  Therefore,  we  have  chosen  (3=1.075  as  the  "universal" 
enhancement  factor;  all  the  calculations  in  [9]  for  the  whole  range  of  c/a  values 
were  done  with  this  p.  No  other  adjustments  of  any  parameters  were  performed. 


\  S  =  2.521 


s  =  2.661 


I  '  '  '  I  ■  ■  ■  I  ■  ■  ■  I  ■  '  ■  I  ■  ■  '  I  ' 


Table  1. 

Dependence  of  the  equilibrium  WS  radius  for  the  FM  BCC  phase  on  the  Stoner 

parameter  enhancement  factor,  p. 


1.000 

1.025 

1.050 

1.075 

1.090 

s,  a.u.a 

2.648 

2.650 

2.654 

2.659 

2.664 

a)  Experimental  value,  s=2.661  a.u.  (Ref.  [19]) 
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We  used  the  LMTO  method  [20]  with  the  so-called  combined  correction 
term  [20],  and  the  Madelung  electrostatic  correction  [21].  Scalar  relativistic 
calculations  on  uniform  meshes  of  1540  points  in  the  irreducible  wedges  of  BCT 
Brillouin  zones  were  done,  with  the  exchange-correlation  functional  of  von  Barth 
and  Hedin  [22].  Also,  the  frozen  core  approximation  was  used  [23]. 

Our  calculations  were  done  for  12  c/a  -values  :  1.0  <  c/a  <  V2.  For  each 
c/a,  self-consistent  non-spin-polarized  calculations  were  performed  for  9  values 

of  the  WS  radius,  s  (  2.521  a.u. 
<  s  <  2.788  a.u.).  In  each 
calculation,  after  convergence 
had  been  achieved,  the  Stoner 
parameters,  /o(c/a,  s),  Eq.(6), 
and  then/  =  1.075/o  were 
found  and  the  averaged  DOS, 
N(m)  was  generated.  Then  the 
Stoner  equation,  Eq.(2),  was 
solved  for  the  equilibrium 
magnetic  moment ,  m,  and  the 
magnetic  energy  Em-  The  data 

base  generated  in  this  series 
of  calculations  will  be  used  in 
our  new  method.  Fig.  3  shows 
the  N(m)  curves  for  BCT 
lattice. 

The  above  data,  however, 
reflect  only  the  tetragonal 
deformation.  An  important 
ingredient  of  the  new  method, 
which  is  to  be  used  in  modeling 
GBs  and  free  surfaces,  should 

0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  be  information  on  free  surfaces. 

The  necessary  data  were 
provided  by  Dr.  R.  Wu  of 
Fig.  3.  The  N(m)  functions  for  FCT  lattices.  Northwestern  University, 

Solid  lines  are  the  BCC  and  FCC  curves.  whose  work  in  Prof.  Art 
The  other  curves  correspond  to  the  Freeman's  group,  has  been  a 

intermediate  c/a  values.  part  of  the  SRG  effort. 

In  the  next  section  we 
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will  describe  in  detail  the  parametrization  of  the  data  and  discuss  the  algorithm  of 
the  method. 


Parametrization  of  the  First-Principles  Data  and  the  Method 
Algorithm. 

Using  the  Stoner  model  with  the  adjusted  exchange  parameter,  we  have 
calculated  the  structural  properties  of  BCC  and  FCC  iron  and,  for  the  first  time, 
the  energetics  of  the  intermediate  states  along  the  Bain  deformation  path  . 

Any  semi-empirical  method  has  to  be  able  to  somehow  describe  the 
environment  of  the  atom  of  interest.  It  is  known  from  a  vast  experience  of  tight- 
binding  calculations  (see,  e.g.[24]),  that  the  number  of  neighbors  is  an  important 
parameter  which  plays  a  crucial  role  in  band-structure  calculations.  We  have 
chosen,  therefore,  to  introduce  the  so-called  "effective  number  of  neighbors", 
Zgff,  to  allow  for  the  atomic  environment.  Since  all  the  quantities  are  to  be 
parametrized  in  terms  of  Zgff,  its  definition  is  not  very  critical.  It  is  important, 
however,  that  Zgff  should  reflect  the  real  atomic  environment.  We  define 

Zeff=2(Rmi„/R)5  exp(l-(R/R„in)2)  (8) 

where  the  summation  is  over  atom  coordinates  R  (Rmin  is  the  distance  from  a 
given  atom  to  its  nearest  neighbor).  No  cut-off  radius  is  introduced,  but,  in  fact, 
the  exponential  provides  a  rather  fast  sum  convergence.  Fig.  4  shows  the 
dependence  of  Zgff  on  c/a  in  body-centered  tetragonal  (BCT)  lattices.  One  can 
see,  that  Zgff  for  BCC  and  FCC  equal  respectively  to  10.37  and  12.63,  while  the 
corresponding  nearest  neighbor  numbers  are  8  and  12.  The  10.37  number 
actually  reflects  the  fact  that  the  6  second  nearest  neighbors  in  the  BCC  lattice  are 
also  situated  rather  close  to  the  first  nearest  neighbors  (R2nd/r^lst=  1-155);  in  the 
FCC  lattice  they  are  farther  away  (R2nd/l^lst=l-414).  Therefore  the  6  second 

nearest  neighbors  in  the  FCC  lattice  are  less  important  giving  a  smaller 
contribution  to  Zgff. 

Apart  from  the  "environment"  or  "structure"  parameter,  Zgff ,  an 
important  parameter  is  the  volume  per  atom.  In  a  perfect  lattice  the  atomic 
volume,  or,  equivalently,  the  WS  sphere  of  radius,  s,  can  be  easily  found.  In  a 
deformed  crystal,  or  a  crystal  with  defects,  the  atomic  volume  is  difficult  to 
define,  and,  in  fact  it  is  a  "bad"  parameter. 
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Fig.  4.  Zgff  (c/a)  for  BCT  lattice. 

As  was  mentioned  above,  the  left-hand  side  of  Eq.  (7),  N(rn),  virtually  does 
not  depend  on  s.  However,  the  right-hand  side  does  depend  on  s,  where  this 
dependence  is  mostly  due  to  the  d-electron  band  width,  W. 

It  is  known  (see,  e.g.,[24]),  that  in  perfect  crystals  W  ~  l/s^.  Therefore,  it 
is  convenient  to  parametrize  W  in  the  form: 

W  =  A(Zeff)(Xl/R'0)i/2  (9) 

As  a  result,  W  has  the  right  volume  dependence. 

The  values  of  W  are  found  as  a  by-product  of  first-principles  calculations: 
our  data  base  contains  the  values  for  a  number  of  BCT  lattices  ,  as  well  as  the 
BCC  (111)  free  surface.  Fig.  5  shows  the  results  of  fitting  the  function,  Eq.(2)  to 
the  calculated  W's.  The  function  A(Zeff)  is  approximated  by  a  polynomial  in  the 
interval  5.2  <  Zgff  <  12.61  (note  the  Zgff  values  between  5.2  and  10.354 
correspond  to  the  (111)  free  surface). 

We  assume  that  the  Stoner  parameter,  I,  is  only  weakly  volume  dependent. 
However,  it  does  depend  on  atomic  environment,  i.e.  on  Zgff.  We  approximated 

the  values  of  /,  Eq.  (6)  (as  found  from  the  first-principles  calculations,  and 
augmented  by  factor  1.075;  see  section  2)  by  two  polynomials  in  Zgff  (Fig.6). 

Thus,  if  N(m)  is  known,  then,  for  a  given  Zgff  and  W,  the  equation  Eq.  (7) 

can  be  solved  for  the  magnetic  moment,  m=ni(W//).  The  next  step  is  calculating 
Ej^in,  the  kinetic  energy  contribution  (the  integral  in  Eq.  (1)),  and  then  the  total 
ferromagnetic  energy,  E^. 
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Y=M0  +  M1*X  +  ...  M8*x®  +  M9*x® 


,  MO 

10.253859995 

Ml 

-0.041673128445 

M2 

0.018964144149 

M3 

-0.0008804504151 

R 

0.99853826057 

10.8 

10.7 

t 

10.6 

10.5 

10.4 


6  7  8  910111213 


Fig.  5.  Factor  A(Zeff)  as  a  polynomial  of  Zgff. 


Y^MO  +  Mrx  +  ...  M8*x^  +  M9*x^ 

MO _ -35.049499939 

M1  23.898350381 

M2 _ -1.6835162297 

M3  _  0.036626169833 

R  1 


x>10.354 


Y=  MO  +  Ml 

*x  +  ...  M8*x®  +  M9*x^ 

MO 

2135.6396608 

Ml 

-528.36686353 

M2 

44.986416556 

M3 

-1.2755515296 

R 

0.9934641922 

Fig.  6.  The  Stoner  parameter,  I,  as  approximated  by  two  polynomials  m  Zgff 


Both  the  A^(m),  m(x),  x=W//,  and  Ekin(m)  can  be  directly  introduced  into 

the  computer  code  as  DATA  statements.  The  three  data  sets  are  the  two- 
dimensional  arrays,  Zgff  being  the  second  variable.  It  is  known,  that 

ferromagnetism  is  unstable  beyond  a  definite  range  of  atomic  volume/structure 
values.  This  translates  into  a  definite  restriction  for  the  values  of  x=W//.  The 
maximal  allowable  value  of  x  as  a  function  of  Zgff  is  shown  in  Fig.7. 

Following  the  above  prescription,  the  computer  code  was  developed,  in  the 
form  of  subroutine  STONFR(ni,  ri,  amag,  emag). 

The  procedure  of  calculations  is  as  follows.  The  subroutine  STONFR  is 
called  from  a  main  program  which  generates  ni  and  ri,  the  total  number  of 
neighboring  atoms  and  their  distances  from  the  atom  of  interest  (ni  should 
include  as  many  atoms  as  it  may  be  necessary  for  the  Zgff  to  converge),  and  amag 

and  emag  are  the  output  magnetic  moment  on  the  atom  of  interest  and  its 
ferromagnetic  energy. 


Fig.  7.  The  maximal  value  of  parameter  x=W/I,  at  which  the  ferromagnetic 
solution  becomes  unstable,  as  a  function  of  Zgff. 

In  the  subroutine,  first  the  Zgff  and  W  (using  function  fac(z)  which  stands 
for  A(Zeff))  are  calculated,  then  the  Stoner  parameter  /(Zgff)  and  W//  are  found 
for  a  given  Zgff.  The  next  step  is  solving  Fq.  (7)  for  the  magnetic  moment  m  - 

simply  reading  the  value  from  the  corresponding  array  using  two-dimentional 
interpolation.  After  finding  m,  the  kinetic  energy,  Fj^jn  (m,Zeff),  is  found  again 
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by  two-dimentional  interpolation  (the  interpolation  subroutine  is  included  in  the 
code).  Fig.  8  shows  the  AV  plots  for  BCC,  c/a=1.24  and  FCC.The  last  step  is 

calculation  of  the  ferromagnetic  energy: 

Em  =  Ekin-l/4/m2  (10) 

Wherever  Eq.  (7)  does  not  have  a  ferromagnetic  solution  (the  ferromagnetic 
solution  is  unstable),  or  is  positive  (the  ferromagnetic  solution  is  metastable), 
the  subroutine  returns  the  values  m=0.,  and  =0. 


Fig.  8.  The  Ekin  AV  plots  for  BCC,  c/a=1.24  and  FCC 


The  test  runs  of  the  subroutine  for  BCT  lattices  for  s=2.661  (equilibrium 
volume  for  BCC  Fe)  from  c/a=1.0  (BCC)  through  c/a=1.32  (at  higher  c/a  at  this 
volume  the  ferromagnetic  state  is  unstable)  showed  an  excellent  agreement 
between. the  values  of  m  and  Em  calculated  by  the  new  method  and  the  results  of 

the  LMTO  calculations  (the  only  exception  is  c/a=1.32,  where  the  energy  is 
already  quite  small).  Table  2  compares  these  results. 


n 


Table  2 

Comparison  of  values  of  m  and  Em  as  calculated  by  the  new  procedure  and  by  the 

LMTO  first-principles  method 


c/a 

Hjcalc 

mLMTO 

(Eb) 

F  calc 

F  LMTO 

(mRy) 

1.00 

2.235 

2.223 

-30.39 

-30.14 

1.08 

2.251 

2.253 

-25.34 

-25.46 

1.12 

2.273 

2.266 

-20.93 

-20.84 

1.16 

2.299 

2.291 

-15.79 

-15.53 

1.20 

2.274 

2.264 

-10.35 

-10.19 

1.24 

2.309 

2.335 

-5.987 

-6.103 

1.28 

2.370 

2.365 

-2.924 

-2.938 

1.32 

2.446 

2.433 

-0.277 

-0.095 

Potential  Applications  of  the  New  Method  in  AFM  and  EAM 

The  new  procedure  can  be  easily  adopted  by  any  semi-empirical  method.  The 
most  popular  method,  the  EAM,  in  its  various  modifications,  fits  some  of  the 
calculated  quantities  to  their  experimental  values.  In  all  EAM  versions,  among  those 
quantities  are  the  bulk  modulus  and  the  three  elastic  moduli;  Cn,  C12  and  C44. 

In  order  to  implement  the  new  procedure,  a  series  of  calculations  for 
hydrostatic,  tetragonal  and  monoclinic  (or  triclinic)  deformations  of  BCC  Fe 
have  to  be  performed,  and  the  ferromagnetic  contributions  to  the  corresponding 
elastic  moduli  found.  Then,  this  information  is  to  be  used  in  fitting  the  adjustable 
parameters  in  the  corresponding  EAM. 

Originally,  it  was  planned  to  implement  the  new  procedure  to  the  recently 
developed  AFM.  The  validity  and  advantages  of  the  AFM  were  then  to  be  tested 
on  the  ferromagnetic  Iron.  However,  the  results  of  our  recent  calculations  [25], 
using  the  AFM  for  modeling  the  relaxation  of  a  grain  boundary  in  tungsten, 
appear  to  be  somewhat  doubtful. 

Our  previous  investigation  using  the  so-called  Finnis-Sinclair  modification 
of  the  EAM  [6]  for  GB  relaxation  in  tungsten,  discovered  an  anomalous 
relaxation  of  the  third  (from  the  GB  plane)  plane  of  W  atoms,  the  so-called  "co¬ 
phase  effect".  It  was  also  observed  in  Fe  by  the  EAM  simulation  [26],  our  Finnis- 
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Sinclair  calculations  [27],  and  restricted  relaxation  calculations  using  the  first- 
principles  method  [7].  However,  this  effect,  which  is  easily  understood  from  the 
physical  point  of  view  [6],  was  not  obtained  in  the  AFM  GB  relaxation  in  W.  This 
negative  result  has  suggested  that  the  detailed  analysis  of  the  AFM  should  be  put 
on  hold.  Instead,  the  emphasis  should  be  given  to  the  implementation  of  our  new 
procedure  to  a  more  traditional  version  of  EAM,  recently  developed  by  the 
author.  This  project  will  be  carried  out  shortly. 
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